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ABSTRACT 

The MEPHISTO space experiments are collaborative United States and 
French investigations aimed at understanding the fundamentals of crystal 
growth. Microgravity experiments were conducted aboard the USMP-1 and 
-2 missions on STS-52 and 62 in October 1 992 and March 1 994 respectively. 
MEPHISTO is a French designed and built Bridgman type furnace which uses 
the Seebeck technique to monitor the solid/liquid interface temperature and 
Peltier pulsing to mark the location and shape of the solid/liquid interface. 

In this paper the Bridgman growth of Sn-Bi and Bi-Sn under terrestrial and 
microgravity conditions is modeled using the finite element code, FIDAP*. 
The numerical model considers fully coupled heat and mass transport, fluid 
motion and solid/liquid phase changes in the crystal growth process. The 
primary goals of this work are: to provide a quantitative study of the thermal 
buoyancy-induced convection in the melt for the two flight experiments; to 
compare the vertical and horizontal growth configurations and systematically 
evaluate the effects of various gravity levels on the solute segregation. Nu- 
merical results of the vertical and horizontal Bridgman growth configurations 
are presented. 


* NASA does not endorse commercial products. Details about the products named in 
this paper were included for completeness and accuracy. No endorsement or criticism of 
these products by NASA should be assumed. 
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1. INTRODUCTION 

The low gravity (g) environment of space has the potential to produce conditions in which 
buoyancy-driven convection is suppressed. This convection-free environment can then be used 
to examine the influence of convection on various phenomena. Of significant interest to industry 
and the scientific community is how convection influences solidification processes. This interest is 
due to the complex relationship among thermal conditions, thermal and solutal convection, growth 
morphology, and macrosegregation in solidified alloy ingots and sensitive electronic materials 
[1,2]. Direct, quantitative analysis of convection in molten metals and semiconductors is particu- 
larly difficult. These materials are generally opaque, which hinders non-intrusive measurements. 
Furthermore, the chemical volatility of the melt makes the insertion of thermocouples impractical. 
Thus, most studies are limited to indirect measurement techniques such as temperature and See- 
beck measurements, Peltier pulsing and post-mortem analysis. The role of numerical modeling 
then becomes crucial in analyzing the system. A better understanding of how convection influ- 
ences solidification can be achieved through synergistic numerical [3] and analytical modeling in 
combination with various solidification experiments on earth [4,5] and in space. 

MEPHISTO [6] is a collaborative effort among Dr. R. Abbaschian (Univ. of Florida at 
Gainesville), Dr. JJ. Favier and teams from Centre d’Etudes Nucleates de Grenoble (CENG) and 
Centre National d’Etudes Spatiales (CNES) in France, and NASA. The MEPHISTO furnace was 
designed and built by the French teams and solidifies three samples simultaneously in a Bridgman 
type arrangement. The French team examined the non-faceted solidification of Sn-Bi alloys in the 
first flight of MEPHISTO. In the second flight of MEPHISTO the faceted growth of the semi-metal 
Bi (alloyed with a small amount of Sn) was studied. The furnace uses the Seebeck technique in 
one of the samples to measure the temperature of the solid/liquid (s/1) interface, Peltier pulsing to 
mark the shape of the interface in the second sample, and has quench capabilities which enable 
the solute profile in the liquid in front of the interface to be determined in the third ingot. During 
flight, all pertinent aspects of MEPHISTO are monitored and controlled through the Payload 
Operations Control Center (POCC), at NASA’s Marshall Space Flight Center, enabling repeated 
melting and solidification runs to be performed at a variety of growth rates. As shown in Figure 
1, the MEPHISTO furnace is located during flight on an experiment bridge known as MPESS 
(the Mission-Particular Equipment Support Structure) in the Shuttle’s cargo bay. Goals of the 
MEPHISTO experiments include determination of how convection influences: a) morphological 
stability of the s/1 interface and the resulting segregation patterns, and b) atomic attachment kinetics 
at the freezing interface, measured in part from the growth rate and undercooling measurements. 

An important part of analyzing processes in space is a full understanding and appreciation for 
the residual gravity present, both quasi steady and transient (p-jitter) [7-10]. For example, during 
MEPHISTO 2 the Shuttle was flown in the -ZLV, +YW attitude, which is with the cargo bay 
toward the Earth and the right (starboard) wing into the velocity vector; and in the -XLV, -ZW 
attitude, which is with the Shuttle’s tail to the Earth and the cargo bay into the velocity vector. These 
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Figure 1. Showing the MPESS carrier, which is located in the Shuttle’s cargo bay, and 
USMP- 2 experiment hardwares: MEPHISTO, Advanced Automated Direc- 
tional Solidification Furnace, the Critical Fluid Light Scattering Experiment 
known as ZENO, Isothermal Dendrite Growth Experiments, and Space Accel- 
eration Measurement System. 

two modes of flight will result in residual gravities of different magnitudes and directions. G-jitter 
is generally assumed to be independent of Orbiter orientation and is measured using SAMS (Space 
Acceleration Measurement System). Steady state residual g can be considered in two parts, i.e. drag 
and tidal acceleration. The g due to aerodynamic drag is always parallel to the velocity of the Shuttle 
and can be estimated from values in the literature which are sensitive to Shuttle orientation (frontal 
surface area) and altitude [7,12]. Tidal accelerations are comprised of rotational or centrifugal 
forces and gravity-gradient contributions [7,9]. Tidal forces vary proportionally with distance from 
the center of mass (CM) of the Shuttle, and are weakly influenced by orbit altitude. Centrifugal 
forces enhance the gravity-gradient acceleration in the local vertical direction, making tidal g 
greater (per unit length) in the direction along the radius of the Earth. The centrifugal force cancels 
the gravity-gradient force in the flight direction. Thus the tidal residual gravity contribution is zero 
in the direction of the Shuttle’s velocity vector, making drag and tidal forces always perpendicular. 
This is significant because, unless an experiment is located at CM, the vector sum of tidal and drag 
forces is always at an angle (other than 0 or 90°) to the velocity vector, making it nearly impossible 
to align the steady state residual g vector with the axis of the experiment sample. In addition, 
the location of an experiment in the cargo bay is not usually known until very late in the flight 
project’s schedule. Thus steady residual g can not usually be considered productively in the design 
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of the experiment hardware. The MEPHISTO furnace on STS-62 was located at approximately 
(x,y,z) = (24,1.05,10.8) meters, with CM at approximately (x,y,z) = (27.7,0,9.4) meters. With these 
locations and literature data [7,12] the tidal forces can be estimated and then summed with the drag 
to determine the total residual steady g and its direction. In this work we examine numerically 
the effects of these quasi-steady accelerations. The detrimental effects of ^-jitter were found to be 
greater than those due to drag and tidal accelerations during MEPHISTO 1. A comparison of the 
quasi-steady and ^-jitter forces is provided elsewhere [12]. Numerical examination of g-jitter and 
Bridgman solidification will be dealt with in a later paper. 

Theoretical work presently being done in support of MEPHISTO includes the examination of 
scaling laws, analytical modeling, and parametric experimental studies [13-23]. Our simulations 
are being used to determine the thermal, solutal and flow fields, and to check the validity, extension 
and generalization of the scaling laws developed. Similar numerical techniques are being used to 
examine other space experiments [23] such as the gallium arsenide crystal growth in microgravity 
experiment [24]. 

In this paper, we report the recent progress in our ongoing numerical modeling in support 
of the MEPHISTO experiments. Our main objectives are: to provide a quantitative study of the 
thermal buoyancy-induced convection in the melt for the two MEPHISTO space experiments; to 
compare the vertical and horizontal growth configurations; and systematically evaluate the effects 
of various gravity levels on the solute segregation. Two basic growth configurations are examined 
in this work: one with the residual gravity vector parallel to the ampoule axis and the second with 
the residual g perpendicular to the ampoule axis. Calculation of the solute field was included in our 
analysis of MEPHISTO 1 but was omitted from our simulations of MEPHISTO 2. And although 
we have solved for interface shape, our presentation of these preliminary results will concentrate 
on the global fluid flow and the likely effects of this flow on macrosegregation. 

2. MEPHISTO GROUND AND SPACE EXPERIMENTS 

All MEPHISTO ground and space experiments have used Bridgman type furnaces with an 
isothermal hot zone, an insulated gradient zone, and an isothermal chill zone. The MEPHISTO 
furnace was designed and built by the French teams and solidifies three samples simultaneously 
in a dual opposed Bridgman type arrangement. Figure 2 shows a schematic of the flight hardware 
and the thermal gradient imposed. The ground (1-g) experiments were done using the engineering 
and flight MEPHISTO furnace duplicates located in France and in similar Bridgman furnaces at 
the University of Florida. All experiments done in France were in the horizontal orientation, with 
both vertical and horizontal Bridgman experiments being performed at the University of Florida. 

All experiments used fused silica ampoules with a 10 mm outer diameter and 6 mm inner 
diameter. Ingots of Sn-0.5at%Bi for MEPHISTO 1 and Bi-0.1at%Sn for MEPHISTO 2 were 
used during flight and in our simulations. The pertinent physical properties for silica, Bi, and 
Sn are tabulated in Table 1 in the Appendix. Some of the values in Table 1 are a result of our 
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Figure 2. Schematic of MEPHISTO furnace and sample temperature profile. 

own distillation of several sources in an effort to obtain more accurate properties over a wider 
temperature range. 

Solidification growth speeds in the ground based experiments ranged from 3.6 to 1 800 rnm/hr, 
and measured thermal gradients at the s/1 interface from 90 to 214 K/cm. During flight, growth 
speeds varied from about 6 to 170 mm/hr with thermal gradients of approximately 200 K/cm in the 
liquid. The thermal gradients imposed along the ampoule wall depend on the temperatures of the hot 
and cold zones and the longitudinal length of the insulated (adiabatic) zone. Information regarding 
the hot and cold zone temperatures and the length of adiabatic zone used in the ground based and 
flight experiments are listed in Table 2 in Appendix. Of critical importance in numerical modeling 
efforts is the accuracy of the input boundary conditions. In this work, the temperatures of the outer 
surface of the ampoule are used as input boundary conditions. However, complete experimentally 
measured temperatures were not available along the length of either the furnace or ampoule. The 
temperatures along the length of the furnace and chill were measured using thermocouples; and, 
ground based MEPHISTO experiments have shown the ampoule surface temperature to be equal to 
the furnace (and chill) temperature at locations away from the gradient region. However, an accurate 
measure of the furnace and ampoule temperature through the adiabatic zone was not available, thus 
the temperature in the gradient zone between the hot and cold zones was calculated using a zero 
radial heat loss condition from the adiabatic zone. Also, in the experiments, the isothermal condition 
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of the hot and cold zones near the adiabatic zone could not be practically maintained due to local 
heat loss at the zone boundaries. To achieve a smoother, more characteristic longitudinal input 
temperature profile, the length of the adiabatic zone was adjusted in the model to give a realistic 
(experimentally verified) temperature gradient in the liquid. 

3. MATHEMATICAL MODEL 

In this paper, we consider heat and mass transport, fluid motion and solid/liquid phase changes 
in the crystal growth process. In particular, the constitutive property of the melt is assumed to be 
Newtonian and its motion is described by the following Navier-Stokes equation: 

Po (^ + u- V«i) = -Vp + V • [p(Vu + (Vu) T )j + A,g[l - A(T - To)] (1) 

where u is the velocity vector, p 0 is fluid density at the reference temperature To, p is pressure, 
fj, is viscosity, T is the temperature variable, g is the acceleration of gravity, (3 t is the volumetric 
expansion coefficient, and the Boussinesq model is adopted to approximate the buoyancy force 
caused by density variation with temperature. Here the change in density associated with species 
concentration has been neglected. The incompressibility condition is given by 

V • u = 0 . (2) 

The heat transport is controlled by the balance of thermal energy 

P0C, (^ + u • V7-) = V ■ (kVT) (3) 

where c p is specific heat and k is thermal conductivity. The solute transport is governed by the 
balance of species concentration 


dC 

dt 


+ u • VC = V • (X?VC) 


(4) 


where C is the concentration of impurity and D denotes the mass diffusion coefficient. On each 
segment of the boundary, it is necessary to prescribe appropriate boundary conditions. Details of 
the computational boundary conditions used in the modeling will be given in section 4. 

The change of phase from liquid to solid is described mathematically by the following phase 
conditions 

Ti(S,t) = T s (S,t) = T m , 
k/VT/ • n - k s VT S • n = p 3 L( u* - u*) • n , 
pi(u l - u*) • n = p s (u s - u*) • n , 
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(u l - u a ) x n = 0 (8) 

which need to be satisfied at the solid/liquid interface. Here subscripts and superscripts, s and l, 
refer to the solid & liquid region, respectively; u a is the solid pulling velocity; u* is the velocity of 
the interface; T m is the melting temperature; n is the unit norm of the interface pointing from the 
liquid to the solid; L is the latent heat of fusion. From a physical point of view, equations (5)-(8) 
represent thermal equilibrium at the interface; the heat flux balance between the phases which 
includes latent heat release; the mass flux balance across the interface and the no-slip condition at 
the liquid side of the interface, respectively. 

In order to include the solute (species) transport in the phase change analysis, it is necessary 
to add two more interface conditions: 


DiVCi • n - D S VC S • n = C s {u s - u*) • n - C t (u l - u*) • n , (9) 

C a = kCi (10) 

where k is the partition coefficient. Eqn. (9) describes the mass conservation for solute transport 
across the interface and eqn. (10) is actually the chemical equilibrium determined by the phase 
diagram. In this case, the melting temperature depends on solute concentration and condition (5) 
becomes 

2) = T s = T m {C) =T a + mCi (5') 

where Ta is a constant and 

m = m(C) = (11) 

is the rate of change of the melting temperature with respect to C . Note that the model summarized 
in this section assumes a sharp solid/liquid interface; the problem will be solved by the interface 
tracking approach with deformable grids. 

4. NUMERICAL SIMULATION 

4.1 The FEM Model 

In this paper we consider two Bridgman growth configurations. The first is the bottom seeded 
(vertical) Bridgman growth. In this configuration the hot (melt) zone is on the top and the cold 
(crystal) zone is at the bottom. The furnace axis is parallel to the gravity vector which is in the 
vertical direction pointing downwards. For this configuration it is reasonable to assume that the 
heat, species and flow fields are all axisymmetric. Therefore a simplified axisymmetric model, 
as depicted in Figure 3, is used to model the vertical Bridgman growth. A 15cm portion of the 
90cm-long rod-shaped sample is modeled. The computational boundary conditions imposed for 
the energy, momentum and species equations are summarized in Figure 4. For the momentum 
equation, we specify no-slip conditions on the surface between the sample and ampoule wall. For 
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heat transport, we impose the measured thermal condition on the outer surface of the ampoule. 
As shown in Fig. 4(b), the thermal profile on the ampoule surface consists of the hot, adiabatic 
and cold zones. Thermal radiation is not included in the model. Our analysis indicates that for 
the material and processing temperature range considered, most of the radiation emitted by the 
heaters gets absorbed by the ampoule and only a small amount of radiation penetrates through the 
ampoule. And, since we impose the furnace temperature on the ampoule wall, radiation is relatively 
less important compared to the longitudinal conduction. The concentration boundary conditions 
are shown in Fig. 4(c). 


Growth Rate V g ► gravity 



Fused Silica Ampoule 

iR 

Z 

Liquid 

Solid 


T 

2 mm 

3mm 




100mm *f* 50mm 


Figure 3. Schematic diagram and geometric definitions of a simplified axisymmet- 
ric FEM model for the vertical Bridgman growth experiment. 


The second configuration is the horizontal Bridgman in which the gravity vector is perpendic- 
ular to the furnace axis (i.e. perpendicular to the 2 -axis defined in Fig. 3). In this case, the solutions 
are assumed to be symmetric about the vertical center plane and an idealized two-dimensional 
model is used to simulate this center plane. Although this 2-D model is clearly a simplification 
of the real situation, we hope it can provide at least some qualitative information for this growth 
configuration. More accurate modeling will rely on the full three-dimensional model which is 
currently under development. 

4.2 Pseudo-Steady-State Model 

The so-called pseudo-steady-state model (PSSM) [25,26] is used to simulate the steady growth 
of BiSn and SnBi in the vertical and horizontal Bridgman configurations. In PSSM, the translation 
of the ampoule is modeled by the melt entering at its hot end with a uniform growth velocity V g 
and composition Co, and by removing the crystal from the cold end at a speed that conserves the 
mass of the alloy in the system. Note that as a simplification of the real growth condition, the 
PSSM neglects the transient effects in the field variables (such as velocity, pressure, temperature 
and concentration, etc.) caused by the steady decrease of the length of the melt during crystal 
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Adiabatic Zone 



dT/dR = 0 

(a) Thermal Boundary Conditions 



u T = 0 u z — free 5 Z 

(b) Velocity Boundary Conditions for Pseudo-Steady State Analysis 



dC/dR = 0 

(c) Boundary Conditions for Solute Concentration 


Figure 4. Computational boundary conditions used for the axisym metric model. 
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growth and the displacement of the ampoule in the furnace. It is known that for sufficiently long 
ampoules the thermal end effects are negligibly small [25]. Another important feature of PSSM is 
the use of the moving coordinate system by fixing its origin at the center of the moving solid/liquid 
interface. 

4.3 The Front Tracking Approach 

The most challenging difficulty posed by phase change problems is the moving solid/liquid 
interface whose position is usually an unknown function of time and space and needs to be 
determined as a part of the solution. In the literature, various numerical techniques have been 
proposed to deal with the moving interface for phase change problems [27]. Among them, two 
main classes of methods can be distinguished, namely the fixed-grid enthalpy methods and the 
front-tracking methods. The front-tracking technique with a deforming mesh is used in our 
modeling. 

Unlike the enthalpy method, the phase change front tracking method can model the phase 
change problems with a sharp (single) solid/liquid interface. The front tracking approach used in 
this work involves a deformable spatial mesh in which nodes located on the interface are allowed 
to move such that they remain on the moving boundary. For each node on the moving interface, 
an additional degree of freedom is introduced. This new degree of freedom determines directly the 
position of the node in space and is an integral part of the representation of the moving interface. 

To update the interface position and remesh the interior domains, a method of spines is used. It 
is a generalization of the method developed by Saito and Scriven [28], in which the moving nodes 
lie on and the interface movement is guided by the generator lines called spines. In particular, the 
position of the moving node is represented parametrically by 

Xi = -I- uti(hj -|- Ac 

Vi = + wti(hj fyj+i)] 4- f3y (12) 

Zi = a z [h j+l + uti(hj - + p 2 

where hj is the interface location parameter for a given spine, (a x ,a y , a z ) is the direction vector, 
(Ac, /3 y , (3 Z ) is the base point of the spine. The location (x*, y,, z z ) of the moving node on the spine 
is determined from its relative position, ujti, to the moving interfaces located at hj and hj + j . Here 
hj are the new degree of freedom introduced. 

For steady state problems, applying Galerkin’s formulation and the standard discretization 
procedure to the momentum equation (1) results in a nonlinear algebraic system in the following 
matrix form [29] ' 

A(U)U + K(T, U) - CP + BX = F (13) 

where X is the global vector of the moving interface unknowns, A(U) accounts for contribution 
from the convective terms, K(U) includes the diffusive terms, C is the matrix discretization of the 
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divergence, B represents the contribution of the normal stress balance condition at boundary, F is 
a body force vector. 

4.4 The Segregated Solver 

System (13) along with the energy equation, species equation and other interface and boundary 
conditions are solved by a segregated solution procedure. In contrast to the traditional fully coupled 
approach (such as the Newton-type solver), the segregated solution algorithm does not form the 
global system matrix directly. Instead, it decomposes the global system matrix into smaller sub- 
matrices each governing the nodal unknowns associated with only one conservation equation. 
These segregated sub-matrices are then solved in a sequential manner. Since the storage of the 
individual sub-matrices is considerably less than that needed to store the global system matrix, 
the storage requirements of the segregated approach are substantially less than those of the fully 
coupled approach. 

4.5 Special Formulation for Concentration 

A special difficulty in solving mass transfer in phase-change problems is the discontinuity of 
solute concentration itself across the interface, as indicated by interface condition (10). Further- 
more, the melting temperature is no longer a constant; instead it has to be determined by the local 
solute concentration. However in the standard FEM formulation, both the trial and test function 
have to be at least C° continuous, and no jump is allowed for the primitive variables. Therefore, a 
special treatment is necessary for the species equation. 

To eliminate the discontinuity of C at the interface, the following transformation 

C* a = C s /k (14) 

is performed for concentration in the solid region. The interface condition (10) is then reduced to 

Ci = C* (15) 

and thus C becomes continuous across the phase change interface. It can be proved [30] that the 
mass flux balance condition (9) can be rewritten in the following form 

(DiVCi - D S VC S ) ■ n = h c C t (16) 


where the coefficient 

h c = p s (l-k)(u s -u*)-n + D s (l-k)n- VC* S /Ci. (17) 

The second term in (17) involves spatial gradient of concentration at the solid side. Since the solute 
diffusivity in solid, i.e. D s , is usually several orders of magnitude smaller than that of the liquid 
Di, it is reasonable to assume this term is negligible. In our computation, the interface condition 
(16) is imposed through a species transfer boundary element at the moving interface. 
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4.6 Numerical Solution 

The axisymmetric and 2-D FEM models are built with the 4-node bilinear element, in which 
velocity, temperature and species are approximated by bilinear shape functions. The pressure is 
approximated as piecewise constant. 

In this work, numerical solutions are obtained by using a modified version of the finite element 
program FTDAP. The details of the FEM formulation in FIDAP are documented in [29], The 
nonlinear iteration termination is controlled by a specified tolerance of 10 -3 for the relative error 
norms of velocity, residual and free surface update. 

5. PRELIMINARY RESULTS FOR MEPHISTO I 

In this section we present our numerical results for the first flight experiment. 

5.1 Two-Step Solution Approach 

In order to simulate the solute segregation, a two-step solution approach was used in our 
modeling. This two-step approach is based on the assumption that the effect of solute field on the 
heat transport and fluid motion is negligibly small. Though the validity of this assumption still needs 
further verification, the main reason for this two-step solution process is the current limitation on the 
front tracking approach in FIDAP. Due to this limitation, it is impossible to include the conjugate 
heat transfer, the moving front tracking and the solute segregation all together in one single model. 

The two-step solution approach can be briefly described as follows. In step one, we consider a 
full ampoule model which includes the solid and liquid SnBi sample as well as the glass ampoule 
wall. A fully coupled heat transport and fluid flow in the liquid SnBi, phase change and the 
conjugate heat transfer in the ampoule wall are solved simultaneously on the basis of this full 
ampoule model. Then in step two, we consider the SnBi sample only (without the ampoule wall). 
On the surface of the SnBi sample, we impose the temperature profile obtained from the first step 
solution. In this way, the effect of the conjugate heat transfer considered in the full ampoule model 
is passed to the simplified model in the second step. The solute segregation coupled with heat 
transport, fluid motion, non-isothermal phase change and mesh deformation are then solved in the 
second step. 

5.2 Buoyancy-Induced Convection 

To study the effects of gravity on the fluid motion in the melt, we plot the contours of stream 
function in Fig. 5 for the horizontal growth at various gravity levels ranging from lO -6 # to 10“ V- 
Note that the plots shown in Fig. 5 are based on the total velocity which includes a constant 
translation velocity and the buoyancy-induced convective velocity. At 10 ~ 6 g, the streamlines 
shown are straight, which indicates that the convective part of the velocity is almost zero. At 10 -5 g 
the streamlines become slightly curved, but the convective velocity is still very small (about 30% 
of the growth rate as indicated by the value of u max ). This suggests that when gravity < 10 -5 g the 
growth is basically diffusion-controlled and the thermal buoyancy-induced convection is negligible. 
At the gravity level of 10~ 4 g, a flow cell has formed in the melt and the flow cell grows bigger 
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(a) 10~ 6 p ti ma x = 1.03901^ 

=T 


(b) 10 -5 <7 u max = 1.3275^ 



(c) 10 - 4 5 u max = 4.3430V,, 



(e) 10- 2 5 u max = 325.75V ff 



(e) 10 -1 5 tt m ax = 3035.51^ 


Figure 5. Contours of stream functions for the horizontal growth with MEPHISTO 
I at various gravity levels. Only the liquid part of SnBi sample is plotted and the 
curved surface on the right is the solid/liquid interface. Here the calculation of 
stream function is based on the total velocity which includes the growth rate V g and 
the buoyancy-induced convective velocity. Namely, u = u c + V g , where u c is the 
convective velocity. 
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with the increase in gravity. The bigger flow cell corresponds to stronger convection. The shallow 
cavity flow pattern shown in Fig. 5 agrees well with the flow pattern reported in the literature for 
the horizontal Bridgman growth [24]. 

The strength of convection in the melt can be quantified by compa rison to t he ma ximum 
velocity. The effect of gravity on the convection in the melt can be seen more clearly from Fig. 6 
in which the total maximum velocity, U max = u^ x +V g , versus the magnitude of gravity is plotted 
for the two growth configurations. The convective strength level marked by the horizontal dashed 
line in Fig. 6 is important, because at this level the buoyancy-driven convective velocity maximum 
is the same as the translation velocity (growth rate). It has been shown that radial segregation is 
often maximum when the convective velocities and the growth rate are similar. Therefore, it can 
be considered as a critical point for the growth of convective strength in the melt. Let g c denote the 
gravity level corresponding to this critical point. For horizontal growth, g c zi 2 x 10 -5 y. When 
gravity is less than g c , the convection is very weak and grows slowly with increases in g. However, 
when the gravity level exceeds g c , convection grows concomitantly with g, and the effects of gravity 
on the fluid motion in the melt become significant. 

For the vertical growth configuration, the buoyancy-induced convection is much weaker than 
for the horizontal case. It is seen from Fig. 6 that for vertical Bridgman, g c « 10 _3 p, and at the 
same gravity level convection is about one to two orders of magnitude lower than that of horizontal 
growth. It is interesting to note that when gravity increases 1000 times from 10~ 3 g to lg, U max 
also grows about 1000 times. The results presented in Fig. 6 indicate that, above the critical value 
g c , the effect of gravity on the convection in the melt is very significant 

5.3 Solute Distribution in the Melt 

For a given bulk alloy, the solute concentration at the interface is affected by the removal rate 
of excess solute atoms rejected at the moving interface. Under ideal conditions (no convection), 
the rejected solute is transported by diffusion only. In reality, however, the transport of excess 
solute atoms to and from the interface is greatly affected by additional factors such as the buoyancy 
force caused by density gradient. In this work, the density variation with temperature is modeled 
by the Boussinesq approximation given in eqn. (1). For SnBi, a positive change of AT will 
reduce the density in the buoyancy term in the momentum equation. The segregation coefficient 
k< 1 and p so iute > Psoivent imply that heaver solute is rejected at the interface. Consequently, the 
double-diffusive convection is solutally stable in a classical vertical growth system. 

For MEPHISTO 1, the initial concentration of Bi in the SnBi alloy was Cq — 0.5 at.%. As 
predicted by the convective strength, the growth is diffusion-dominated at the microgravity level 
of 10~ 6 p. To verify the numerical results, we compare the simulated solute distribution on the 
centerline of the 2-D horizontal growth model against the available 1-D analytical solution. It 
is seen from Fig. 7 that the modeled solute boundary layer at 10 _6 p is indeed very close to the 
diffusion-controlled solution. 

In order to describe quantitatively the transition from diffusion-controlled growth (without 
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Figure 6. Effects of gravity on the strength of convection in the melt 
measured by the maximum total velocity, C/ max > which includes the growth 
rate V g and the thermal buoyancy-driven convective velocity. 



Figure 7. Comparison of the simulated axial solute distribution on the cen- 
terline of the SnBi sample with the 1-D diffusion-controlled growth solution. 
The FEM result is based on the 2-D model for the horizontal configuration. In 
order to display the thin boundary layer, only a portion of the entire centerline 
is shown. 
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bulk convection) to growth with intensive laminar convective mixing, we compute two -important 
segregation parameters. The first is the effective segregation coefficient fc e ff defined as [26] 

fc eff = k <C>i / «C» (18) 

where < C >/ is the lateral average of solute concentration over the solid/liquid interface and 
C » is the volumetric average of solute concentration in the melt The second is the percentage 
radial segregation AC defined as 

AC = |£C| max /C 0 (19) 

or 

AC = |C top — Cbottoml/Co . (20) 

Here |5C| max is the maximum difference in liquid concentration at the interface. For directional 
solidification, the diffusion-controlled growth with a planar interface leads to uniform radial solute 
distribution, thus we may have AC « 0. However, if the s/1 interface is not flat for example, 
due to radial temperature gradients, significant radial solute gradients may develop even under 
diffusion-controlled conditions. If the melt is sufficiently long that the diffusion layer near the 
interface occupies only a small fraction of the total length, fc e ff approaches unity when there is no 
convective flow in the melt other than the unidirectional growth velocity V g . The other limit of k e fi 
represents the steady-state well-mixed growth in which the intense convection leads to a complete 
mixing in the melt. In this case the value of k e ff approaches k. 

The computed fc e ff and AC for both the vertical and horizontal growth are plotted as a function 
of gravity levels in Figs. 8 and 9, respectively. For the vertical growth case, the values of k e ff are 
very close to one when gravity is less than 10 -2 <? as shown in Fig. 8(a). The value of k e ^ starts to 
decrease rapidly at 10 ~ 2 g and does approach k with increasing gravity levels. The corresponding 
radial segregation presented in Fig. 9(a) shows that there is a radial solute gradient at the interface 
even in the diffusion dominated growth range. The radial segregation reaches a maximum at about 
10 _2 <? and then decreases with the increase of gravity. The radial segregation in the liquid is caused 
by the nonplanar shape of the interface, and since we feel improvements are still possible in our 
interface shape results, we will be examining segregation in more detail in a later paper. For the 
horizontal configuration, the value of k e ff starts to decrease dramatically at a much lower gravity 
level of 10~ 4 <7, because of the much stronger convection. It is seen from Fig. 8(b) that k e ff is 
very close to k at 10 _1 <? which indicates that a well-mixed state has been achieved at this gravity 
level. Note that although the values of the radial segregation computed by eqns. (19) and (20) 
are different, the predicted maximum AC occurs at the same gravity level of 10 -4 <? for horizontal 
growth as shown in Fig. 9(b). The plots of fc e ff and AC vs. the magnitude of gravity presented 
here agree very well with the schematic curves given by Kim, Adomato and Brown in [26]. 

By considering the worst orientation of gravity vector during the space experiment (namely the 
horizontal configuration), the results presented in Figs. 8 and 9 suggest that the growth is basically 
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Gravity 


Figure 8. Effects of gravity on the axial solute segregation measured in 
terms of the effective segregation coefficient, fc e ff- (a) For the vertical growth; 
(b) for the horizontal configuration. 
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diffusion-controlled when the magnitude of gravity is less than lO -5 ^. Therefore our numerical 
observation supports the diffusion-controlled growth assumption used by the French team [13,16]. 
Another interesting point is the gravity level at which the radial segregation reaches its maximum. 
For both growth configurations, the maximum AC occurs at the approximate ending point of the 
gravity range over which k e fi indicates diffusion-controlled growth. For example, for the horizontal 
growth the ending point of the gravity range for fc e ff ~ 1 is 10 ~ 4 < 7 , that is exactly the same gravity 
level at which AC has its maximum, as one can see from Fig. 8(b) and Fig. 9(b). Actually this 
gravity level corresponds to a weak convection in the melt with an incomplete mixing of solute. 
This in turn leads the concentration field adjacent to the interface highly distorted and the amount of 
radial segregation reaches its maximum, though the axial segregation (measured by fc e ff) indicates 
diffusion-controlled growth. 

5.4 Convergence of Numerical Solutions 

In order to check the convergence and accuracy of the numerical solutions, a set of three 
meshes were used in the modeling for the vertical growth case. The three meshes are referred to as 
mesh 1 to 3 in Fig. 8(a) and Fig. 9(a). Mesh 1 is the coarsest with a total of 820 bilinear elements 
and 913 nodes. Mesh 2 is the intermediate with a total of 1344 elements and 1469 nodes. Mesh 3 
is the finest with about 3000 elements. 

For velocity, the numerical solutions from the three meshes are all very close. For example the 
difference between the maximum velocities predicted by mesh 1 and mesh 3 is less than 1 %, which 
suggests that even mesh 1 can provide a very good solution for velocity. Our analysis shows that 
the most difficult variable for numerical convergence is the solute concentration, perhaps due to 
the very thin solute boundary layer near the interface. For the effective segregation coefficient, the 
results presented in Fig. 8(a) show that when gravity is less than 10 -1 <jr the values of k e fj computed 
from the three meshes are very close. But the error of the mesh 1 becomes unacceptable at lg. 
For the calculation of the radial segregation, mesh 1 fails completely as can be seen in Fig. 9(a). 
This suggests that a coarse mesh is not capable of resolving the thin solute boundary layer. In all 
the cases, the differences between the solutions obtained from meshes 2 and 3 are very small. Our 
experience also indicates that grading towards the interface (i.e. gradually reducing element size 
towards interface) helps improve the solution accuracy for solute concentration. 

6. PRELIMINARY RESULTS FOR MEPfflSTO H 

For MEPHISTO 2, seven gravity levels were considered: lO^g, 10' 5 g, l(T*g, 10'^g, 10 2 g, 
10‘ 2 g and lg. The 10' 6 g to 10' 3 g cases simulate the range of accelerations expected during space 
processing and the lg case simulates ground based experiments. 

Figs. 10(a), 10(b) and 10(c) show the velocity vectors for three gravity levels for the vertical 
configuration. It is seen that the flow cell forms at a g-level of 10' 5 g and builds up as gravity in- 
creases. At a gravity level of 10' 3 g, the flow velocities become comparable to the growth velocity. 
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Figure 10. BiSn, vertical configuration: convective cells at (a) lO^g with u c TOLC =7.53xlO' 7 
cm/sec, (b) lCMg with u c TOK =7.72xlO' 4 cm/sec, and (c) lg with \i c max ^2.59xl0' 1 cm/sec. 


20 





M. Yao, R. Raman and H.C. de Groh HI 



Figure 11. BiSn, Horizontal configuration: Convective cells for (a) 10°g with 
u‘ =2.02xl0- 5 cm/sec, (b) 10' 5 g with u c max =2.04xl0‘ 2 cm/sec, (c) lg with u c max =3.72 . 
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A very strong convective cell is seen for the lg case. Figs. 1 1(a), 1 1(b) and 1 1(c) show the veloc- 
ity vectors for three gravity levels for the horizontal configuration. It is seen that the velocities are 
greater and stronger convective cells are formed for this configuration. 

To clearly correlate the effect of gravity on the flow velocities, the ratio of the maximum ve- 
locity to the growth velocity ~V g was plotted against each g-level in Fig. 12. This plot shows that 
the log of the velocity ratio is nearly linear with the log of the gravity. This means that 
grows with the same order of magnitude as the gravity level. Fig. 12 also compares the growth of 
uC max f° r the horizontal configuration. It is interesting to see that the slope is nearly the same for 
this case, showing that the flow velocities are consistently an order of magnitude higher for the 
horizontal configuration. 

It has been shown here and by others 2 that the greatest amount of radial segregation occurs 
when convective velocities near the interface are close to the solidification rate (u c m<r yv^ approx- 
imately = 1). Fig. 12 shows that for the vertical case, convective levels are expected to be most 
damaging to the radial homogeneity at about 10'^g, because convective velocities are about the 
same as the growth velocity. 

For the horizontal case however, concern for radial homogeneity is justified at about 3x1 O' 5 
g. The steady state background g-levels on the Space Shuttle are generally accepted to be less 
than 10" 5 g, thus steady g-levels imply diffusion controlled growth for MEPHISTO-2 since con- 
vective levels are expected to be significantly less than the growth velocity. 



Figure 12 . A log-log plot of u c m ^/V^ vs. g-levels showing the effect of grav- 
ity on convective maximum velocities for MEPHISTO 2. 
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Residual gravity transients caused by thruster firings and other disturbances are often in the 
range of lO^g and greater for short periods; the response of solidification processes to non-steady 
gravity is presently under investigation and will be dealt with in a later paper. However, results 
from the first flight of MEPHISTO show experimental indications of convection during, for ex- 
ample, crew exercise, when the magnitude of the transient acceleration is in the range of 5x1 O^g. 

The sample temperature in the hot and cold zones and the profile through the adiabatic zone, 
which was calculated as part of the solution, are shown in Fig 13. The temperature contours for 
two gravity levels for the horizontal configuration are shown in Fig 14. It is clear that the contours 
become distorted by the flow field, with warmer liquid traveling along the top wall and cooler liq- 
uid returning along the bottom. For the vertical configuration, the distortion is much less because 
the velocities are lower. However, since at 1-g in the horizontal orientation, the interface is ex- 
pected to be curved (with the solid concave, or female) 5 further modeling of the interface shape is 
presently underway. 

During the first flight of MEPHISTO, under some conditions, convection in the melt was in- 
dicated by fluctuations in the interface temperature (as measured using the Seebeck technique). 
By nondimensionalizing these conditions we can project convective conditions to MEPHISTO 2, 
and in combination with the numerical modeling results, estimate the relative influence of con- 
vection and diffusion during MEPHISTO. The Prandtl number, Pr = ~ , compares the rate of 
diffusion of momentum and heat where v is the dynamic viscosity and a is the thermal diffusivity. 
The Schmidt number. Sc = ^ , similarly compares diffusivities of momentum and species where 


xlO 5 



Figure 13. BiSn, calculated temperature profile at the center-line of the sample 
(R=0) for the vertical configuration based on the 2D axisymmetric model. 
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Key: 

A: 326.5 K 
J: 749.5 K 
AT = 47 K 



(b) GRAVITY = lg 

Figure 14. MEPHISTO 2 temperature contour plots for (a) lO^g and (b) lg in the hor- 
izontal configuration. The temperature difference between isotherms is 47K. 


D is the species diffusivity. The thermal Grashof number is given by: 

_gP 1 ATL i ' 

Gr T - 2 

v 

where is the coefficient of thermal expansion, AT = T h - T m is the temperature difference be- 
tween the hot zone, T h , and the alloy melting temperature, T m . The distance over which AT is ac- 
tive is L, the characteristic length of the system, i.e., the longitudinal distance between T h and T m 
(L»3.7cm for MEPHISTO I, L»2 cm for MEPHISTO 2). The solutal Grashof number is given by, 

gP c AC5 3 

0r * = — 

V 

where |3 C is the solutal expansivity, the characteristic length, 8, is the boundary-layer thickness 
and AC is the concentration difference between the far field liquid and the liquid at the interface. 
The solutal Peclet number is given by Pe = u c max 5/D, where u c max is the maximum convective 
flow velocity in the melt Table 1 shows the dimensionless groups for MEPHISTO 1 and 2. 
Prandtl and Schmidt numbers were evaluated at 20 K above T m . 

The thermal Grashof numbers for the MEPHISTO experiments are of the same order, thus 
convective mixing is expected during crew exercise (at g levels above about 300 jig) and during 
such activities as major truster firings (at 2500|iig) as was observed during MEPHISTO 1. 
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Amounts of convection insufficient to cause disturbances detected by Seebeck interface tempera- 
ture measurements may still influence the solidification process significantly. For example, rather 
slow flow velocities have been shown to be most detrimental to radial homogeneity (Brown 
(1988), Yao and de Groh (1994)). The solutal Peclet numbers shown in Table 1 indicate that con- 
vective and diffusive transport are of the same order when gravity is in the 10 to 100 pg range. 
This indicates that for MEPHISTO 1 and 2 diffusion-dominated growth was likely during sleep 
periods (at g levels less than about 5pg). This is consistent with Fig. 9 results, which show the ra- 
dial solutal gradient in the liquid to be relatively small at g levels below 10 ' 5 g. 


7. CONCLUSIONS 

The goal of this work was to, using numerical techniques, quantify buoyancy-driven convection 
during Bridgman solidification and to determine the effect of gravity on convective velocities in the 
vertical and horizontal configurations for two recent NASA space experiment projects, MEPHISTO 
1 and 2. For MEPHISTO 1, the solute field was included in our analysis, however, for MEPHISTO 
2 , solute was omitted. ... 

Numerical simulation indicates that the buoyancy-induced convection in the melt intensifies 
quite rapidly as gravity increases. At equivalent gravity, flow velocities occurring during horizontal 
conditions were about 15 times greater than those during vertical growth for both MEPHISTO 1 
and 2. By considering the worst orientation of gravity vector during the space experiment (namely 
the horizontal configuration), our results suggest that the growth is basically diffusion-controlled at 
^-levels below 1 O “ 5 0 for both MEPHISTO 1 and 2. Therefore our numerical observation supports 
the diffusion-controlled growth assumption [13] and the assumption that under quasi-diffusive 
conditions residual convection does not perturb the concentration field significantly at low 0 -levels 
[16] used in the analysis by the French team. 
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Table lPhysical Property Data for amorphous Si0 2 ampoule, Bi, and Sn. 31 ' 34 

Property 

Ampoule Bi solid 

Bi liquid Sn solid 

Sn liquid 

Density, p, g/cc 

2.20 9.80 

10.07 7.30 

7.00 

Melting point, T m , °C 

271.3 

231.912 


Heat of fusion, AH, cal/g 

12.5 

14.5 

£ 



Thermal expansion coef., |3, K 1 -125x10 6 -88xl0' 6 

Diffusivity of Sn in liquid Bi: D Sn ^ B i = 5.2xlO' 4 (exp(-3.2kcal-mol' 1 /RT) cm 2 /s 
Diffusivity of Bi in liquid Sn: D Bi inSn = 3.355xl0- 4 (exp(-2.98kcal-mol 1 /RT) cm 2 /s 
Viscosity, ji,, cP: Liquid Bi: Ji. = 0.4458cP exp[ 1.541 kcal-mol 1 /RT] 

Liquid Sn: ji = 0.5382cP exp[1.3 kcal-mol J /RT] 

Partition coefficient for Bi rich alloy with Sn, k = 3.5xl0' 3 at%/at% 

Partition coefficient for Si rich alloy with Bi, k = 0.29 at%/at% 

Thermal conductivity, k, W/cm-K: 


at 300 K 


0.0786 


0.666 


400 K 


0.0704 


0.622 


423 K 

0.0189 





473 K 

0.0192 





500 K 


0.0663 




505 K 




0.595 


506 K 





0.303 

544 K 


0.065 




545 K 



0.124 



573 K 

0.0201 





600 K 



0.131 


0.323 

673 K 

0.0213 





700 K 



0.141 


0.343 

773 K 

0.0222 





800 K 



0.15 


0.364 

873 K 

0.023 





900 K 



0.159 


0.384 

973 K 

0.0239 





Specific heat, CL, cal/g-K: 






at 300 K 

0.17 

0.0295 


0.055 


400 K 

0.21 

0.031 


0.058 


500 K 

0.235 

0.033 


0.063 


506 K 





0.064 

520 K 


0.034 




534 K 


0.043 




538 K 


0.0567 




545 K 



0.0346 



600 K 

0.25 


0.0336 


0.0586 

700 K 

0.26 


0.0326 


0.0568 

800 K 

0.27 


0.0321 


0.055 

900 K 

0.28 
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Table 2 Approximate furnace hot and cold zone temperature ranges in °C, adiabatic zone lengths 
in mm, angle off vertical between g and the ampoule axis in degrees - for experiments in space the 
angle is approximate and variable due to changes in flight mode, and magnitude of gravity level, 
on Earth and in space, as well as the particular growth speeds simulated, in units of mm/hr. 

Experiments on: Tw T^y Adiabatic width Growth speed 

Angle 

g 

Earth 500 to 760 30 

In Space 600 to 760 50 

9 to 14 
14 

6 to 1800 
6.67 to 144 

0 and 90 
50 to 70 

1 

lO^to 10’ 3 

MF.PWTSTO 1 simulations: — 

500 50 (30@lg) 

20 

6.67 

0 and 90 

lO^to 1 

MPPTTTSTO 5 simulations: — — — 

Vertical 500 30 

18 

20.5 

0 

lO^to 1 

Horizontal 500 30 

20 

20.5 

90 

lO^to 1 


Table 3 Dimensionless numbers for MEPHISTO 1 and 2. For MEPHISTO 1, convection was 
observed at lg, 2500 pg, and 300 pg. When the signs of the Gr s and Gr T are the same the solutal 
and thermal effects tend to enhance each other. U p is growth velocity in units of mm/hr. 

Dimensionless 

MEPHISTO 1 

MEPHISTO 2 

Number 

Up = 7 Up = 100 

Up = 7 Up = 100 

Pr 

0.016 

0.018 

Sc 

140 

60 

Gr s 

1 pg 

-1.3X10- 4 -4 x 10" 8 

0.09 2.8xl0‘ 5 

300 pg 

-0.04 -1.2xl0* 5 

27 8.4xl0' 3 

2500 pg 

-0.33 -l.OxlO- 4 

226 0.07 

1 g 

-130 -0.04 

9-OxlO 4 28 

Gr r 



1 pg 

-23 

-32 

300 pg 

-7.x 10 3 

-9.7X10 3 

2500 pg 

-5.7xl0 4 

-8.1X10 4 

lg 

-2.3x1 0 7 

-3.2xl0 7 

Pe @ lOOpg 

7 

20 


30 
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